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Abstract. Spatiotemporal correlations of the one-dimensional spring-block (Burridge- 
KnopofF) model of earthquakes, either with or without the viscosity term, are studied 
by means of numerical computer simulations. The continuum limit of the model is ex- 
amined by systematically investigating the model properties with varying the block-size 
parameter a toward a — > 0. The Kelvin viscosity term is introduced so that the model 
dynamics possesses a sensible continuum limit. In the presence of the viscosity term, many 
of the properties of the original discrete BK model are kept qualitatively unchanged even 
in the continuum limit, although the size of minimum earthquake gets smaller as a gets 
smaller. One notable exception is the existence/non-existence of the doughnut-like qui- 
escence prior to the mainshock. Although large events of the original discrete BK model 
accompany seismic acceleration together with a doughnut-like quiescence just before the 
mainshock, the spatial range of the doughnut-like quiescence becomes narrower as a gets 
smaller, and in the continuum limit, the doughnut-like quiescence might vanish altogether. 
The doughnut-like quiescence observed in the discrete BK model is then a phenomenon 
closely related to the short-length cut-off scale of the model. 



1. Introduction 

An earthquake is a stick-slip dynamical instability of a 
pre-existing fault driven by the motion of a tectonic plate 
[Scholz, 1998, 2002]. While an earthquake is a complex 
phenomenon, certain empirical laws such as the Gutenberg- 
Richter (GR) law and the Omori law concerning its statisti- 
cal properties are known to hold. One fruitful strategy in un- 
derstanding such statistical properties of earthquakes might 
be to properly model an earthquake fault and to elucidate 
the model properties by numerical computer simulations. 
One of popular statistical models of an earthquake fault 
is the so-called spring-block model originally proposed by 
Burridge and Knopoff (BK) [Burridge and Knopoff, 1967]. 
In this model, an earthquake fault is simulated by an as- 
sembly of blocks, each of which is connected via the elastic 
springs to the neighboring blocks and to the moving plate. 
All blocks are subject to the friction force, the source of 
the nonlinearity in the model, which eventually realizes an 
earthquake-like frictional instability. While the spring-block 
model is obviously a crude model to represent a real earth- 
quake fault, its simplicity enables one to study the statistical 
properties with high precision. 

Carlson, Langer and others [Carlson and Langer, 1989a; 
Carlson and Langer, 1989b; Carlson et al., 1991; Carlson, 
1991a; Carlson, 1991b; Shaw et al, 1992; Carlson et al., 
1994; Schmittbuhl et al, 1996] studied the statistical prop- 
erties of the BK model quite extensively both in one dimen- 
sion (ID) and two dimensions (2D), paying particular at- 
tention to the magnitude distribution of earthquake events. 
In these studies, a simple velocity-weakening friction law 
was used as a constitutive law. The spring-block model has 
also been extended in several ways, e.g., taking account of 
the effect of the viscosity [Myers and Langer, 1993; Shaw, 
1994; De and Ananthakrisna, 2004], modifying the form of 
the friction force [Myers and Langer, 1993; Shaw, 1995; De 
and Ananthakrisna, 2004], considering the long-range inter- 
actions between blocks [Xia et al., 2005, 2008a; Mori and 
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Kawamura, 2008b], driving the system only at one end of 
the system [Vieira, 1992], or by incorporating the rate- and 
state-dependent friction law [Ohmura and Kawamura, 2007]. 

In the previous papers, the present authors also per- 
formed an extensive numerical study of the statistical prop- 
erties of the ID and 2D BK models with the velocity- 
weakening friction law for a wide parameter range, focusing 
on their spatiotemporal correlations [Mon and Kawamura, 
2005; 2006; 2008a; Kawamura, 2006]. These studies have 
revealed several interesting features of the BK model. For 
example, the model exhibits either "near-critical", "super- 
critical" or "subcritical" behaviors primarily depending on 
its frictional parameter representing the extent of the ve- 
locity weakening. Preceding the mainshock, the frequency 
of smaller events is gradually enhanced, whereas, just be- 
fore the mainshock, it is suppressed in a close vicinity of 
the epicenter of the upcoming event (the Mogi doughnut). 
The apparent _B-value of the magnitude distribution changes 
significantly preceding the mainshock. 

Although the BK model has widely been used as a use- 
ful tool to investigate statistical properties of earthquakes, 
the block discretization inherent to the model construction 
is a crude approximation of the originally continuum earth- 
quake fault. It introduces the short-length cut-off scale into 
the problem. Therefore, in order to check the validity of 
the model, it is crucially important to examine the contin- 
uum limit of the BK model. Indeed, Rice criticized that the 
discrete BK model with the velocity-weakening friction law 
was "intrinsically discrete" , lacking in a well-defined contin- 
uum limit [Rice, 1993]. Rice argued that the spatiotemporal 
complexity observed in the discrete BK model was due to 
the inherent discreteness of the model, which should disap- 
pear in continuum. This problem of the continuum limit 
of the BK model was also discussed by Myers and Langer 
[Myers and Langer, 1993], who introduced the Kelvin vis- 
cosity term to produce a small length scale which allows a 
well-defined continuum limit. Myers and Langer, and subse- 
quently Shaw [Shaw, 1994] , observed that the added viscos- 
ity term smoothed the rupture dynamics, apparently giving 
rise to the continuum limit accompanied by the spatiotem- 
poral complexity. 

The BK model represents an earthquake fault by an as- 
sembly of blocks and springs. Conversely, one can regard 
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the BK model as obtained by discretizing an appropriate 
wave equation describing the continuum. In one spatial di- 
mension, an appropriate wave equation might be given by 
the partial differential equation, 
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where U{x',t') is the displacement at the spatial position x' 
and time i', u' the loading rate by the plate motion, ^' the 
wave velocity, uj the characteristic frequency and <j}' is the 
friction force per unit mass. One may discretize Eq.(l) as 
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where a' is the grid spacing corresponding to the block size 
of the BK model. If one regards Ui as the displacement 
of the i-th block, Eq.(2) is the equation of motion of the 
ID BK model. Conversely, if one takes the limit a' -^ 
starting from the discrete BK model (2) , one can obtain the 
continuum limit of the BK model. 

As mentioned, the naive continuum limit of the discrete 
BK model defined by Eq.(2) has a problem in that the pulse 
of slip tends to become increasingly narrow in width in the 
limit a' —> 0, i.e., the dynamics becomes sensitive to the 
grid spacing a' . One way to circumvent this problem is to 
introduce the viscosity term rj'd^U/dx' dt' into Eq.(2) to 
produce a small length scale. Myers and Danger showed 
that, owing to the added viscosity term, the system became 
independent of the grid spacing a' as long as a new small 
length scale e', defined by 
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is sufficiently larger than the grid spacing a' [Myers and 
hanger, 1993]. Here, r]' is a viscosity coefficient, and a is 
the velocity-weakening parameter representing the rate of 
the friction force getting weaker on increasing the sliding ve- 
locity (to be defined below in Eq.(7)). Shaw also showed by 
adding the viscosity term to the ID BK model that the mag- 
nitude distribution became independent of the grid spacing 
a' for sufficiently small a' [Shaw, 1994]. 

In the present paper, following Shaw, we numerically ex- 
amine the continuum limit of the ID BK model by perform- 
ing extensive numerical simulations on the BK model with 
successively smaller grid spacings a' . Namely, with varying 
the grid spacing a' , we study how seismic spatiotemporal 
correlations of the ID BK model vary and approach the 
continuum limit. The two types of the BK model, i.e., the 
one without the viscosity term (the original BK model) and 
the one with the viscosity term (the modified BK model) 
are studied. Our work can be regarded as an extension of 
the work of Shaw concerning the magnitude distribution of 
the ID BK model to various types of seismic spatiotemporal 
correlations of the same model. 

The present paper is organized as follows. In §2, we in- 
troduce the model and explain some of the details of our nu- 
merical simulation. In §3, we examine the continuum limit 
of the ID BK model without the viscosity term. We analyze 
the magnitude distribution and various types of seismic spa- 
tiotemporal correlation functions, including the recurrence- 
time distribution, the time- correlation function of seismic 
events before and after the mainshock, the time develop- 
ment of the spatial correlation function of seismic events 
before the mainshock, and the time development of the mag- 
nitude distribution function before and after the mainshock. 
In §4, we make a similar analysis for the ID BK model with 
the viscosity term. Finally, §5 is devoted to summary and 
discussion. 



2. The model and the method 

Let us begin with the ID wave equation in the continuum 
with the Kelvin viscosity term rj'&^U /dx' dt' , 
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In order to make the equation dimensionless, we measure 
the time t' in units of the characteristic frequency ui, the co- 
ordinate x' in units of £,' /co and the displacement Ui in units 
of the length £ — (j)'{Q)/uP , 4>'{0) being a static friction per 
unit mass. Then, the equation of motion can be written in 
the dimensionless form as 
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where t = t'u) is the dimensionless time, x — x' /{(,' /uj) the 
dimensionless coordinate, u = U/£ the dimensionless dis- 
placement, 0(m) = (f)' {tf) / {Cup) the dimensionless friction 
force, and r; = r?'/(C' 1^) is the dimensionless viscosity co- 
efficient. 

After the block discretization, one gets 



ut — Ui + -^(ui+i — 2ui + -Ui-i) 



(6) 



where a = a' /{(,' /uj) is the dimensionless grid spacing, and 
the notation u = du/dt has been used. If we put ?) = 
in Eq.(6), it reduces to the BK model used in the previous 
papers with the correspondence I — 1/a [Mori and Kawa- 
mura, 2005; 2006]. Note that the dependence on ^' has been 
adsorbed into the definition of the dimensionless variables. 
In the following, we study the statistical properties of the 
model (6) with varying the dimensionless grid spacing a in 
order to see how the spatiotemporal correlation functions 
behave in the continuum limit a — > 0. We investigate both 
cases of the non- viscosity model with rj = and the viscosity 
model with rj = 0.02. 

In order for the model to exhibit a dynamical instabil- 
ity corresponding to an earthquake, it is essential that the 
friction force possesses a frictional weakening property, 
i.e., the friction should become weaker as the block slides. 
Here, we assume the velocity-weakening friction force. The 
velocity- weakening friction force used by Carlson [Carlson 
et at, 1991] was defined by 
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where its maximum value corresponding to the static fric- 
tion has been normalized to unity. As noted above, this 
normalization condition (j){u = 0) = 1 has been utilized to 
set the length unit £. The back-slip is inhibited by imposing 
an infinitely large friction for iii < 0, i.e., (f>{u < 0) = -co. 
The friction force is characterized by the two parameters, a 
and a. The former, a, introduced by [Carlson et al., 1991] 
as a technical device, represents an instantaneous drop of 
the friction force at the onset of the slip, while the latter, 
a, represents the rate of the friction force getting weaker on 
increasing the sliding velocity. In the present simulation, we 
regard a to be small and fix a = 0.01. 

Previous studies on the model has revealed that the 
velocity-weakening parameter a is most important among 
various parameters of the model in determining its statistical 
properties [Mori and Kawamura, 2005; 2006; 2008a; 2008b]. 
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In the present study, we set a either a = 1 or a = 3 as two 
typical cases. The former a = 1 corresponds to the "near- 
critical" case in which the magnitude distribution exhibits a 
near straight-line behavior close to the GR law, whereas the 
latter case a = 3 corresponds to the "supercritical" case in 
which the magnitude distribution exhibits a characteristic 
peak at a larger magnitude with the GR-law-like behavior 
realized only at smaller magnitudes. 

We solve the equation of motion (6) by using the Runge- 
Kutta method of the fourth order. The width of the time 
discretization At is taken to be At = 10^"^. Total number of 
10® ~ 10^ events are generated in each run, which are used 
to perform various averagings. In calculating the observ- 
ables, initial 10* ~ 10^ events are discarded as transients. 
In order to eliminate the possible finite-size effects, the total 
number of blocks A'' arc taken to be large. A'' = L/a with 
fixing L — 200, periodic boundary condition being applied. 
In a very special occasion where the finite-size effect is most 
severe, we simulate even larger system of L = 400. The 
continuum limit of the model is investigated by varying the 
dimensionless grid spacing a in the range 1 > a > 1/32. The 
total CPU time is about 3,600 hours of of Pentium D. 

The dimensionless distance r between the block i and i' 
is measured by 

r^a\i-i'\. (8) 

The small length scale e' defined by Eq.(3) is also given in 
the dimensionless form as 
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The magnitude of a seismic event, /i, is defined by a loga- 
rithm of its moment M, i.e.. 



^ = In M = In 
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where Au; is the total displacement of the i-th block during 
a given event and the sum is taken over all blocks involved 
in the event [Carlson et al., 1991]. 



magnitude distribution is considerably widened as the grid 
spacing a becomes smaller. The maximum magnitude gets 
bigger and and the minimum magnitude gets smaller. The 
observation that the size of minimum event gets smaller for 
smaller a can naturally be understood, since a sets a small 
cut-off length scale. As can be clearly seen from the figures, 
R[^) does not seem to converge to an asymptotic form even 
at our smallest grid spacing a — 1/32. Namely, R{i-i) contin- 
ues to change its form as a gets smaller, although the change 
seems to be limited to that of the range and the extent of 
R{li) for a < 1/8. Similar results were reported by [Shaw, 
1994]. 

3.2. The local recurrence-time distribution 

Next, we examine the continuum limit of the recurrence- 
time distribution of the model. Following our previous stud- 
ies on the discrete BK model [Mori and Kawamura, 2005; 
2006], we study here the local recurrence-time distribution 
function P{T). The local recurrence time T is defined by 
the time until the next event occurs with its epicenter lying 
in a vicinity of the previous event within the distance r from 
the epicenter of the previous event. We consider events of 
their magnitude fi > fic, and compute the distribution of 
the local recurrence time T with r — 7.5. The magnitude 
threshold is set to /ic = 2. 

In the main panel of Fig. 2, we show on a log-log plot the 
computed local recurrence-time distribution function P{T) 
for the cases of a = 1 (a) and oi^a = 3 (b). The recurrence 
time is normalized by its mean T, which is Tv — 2.22, 0.56, 
2.30, 2.67, 2.59 and 3.20 (respectively for a = 1, 1/4, 1/8, 
1/12, 1/16 and 1/24) for a = 1 (Fig. 2(a)), and Tu = 41.4, 
1.72, 2.76, 3.90, 4.81, 6.09 and 6.72 (respectively for a = 1, 
1/4, 1/8, 1/12, 1/16, 1/24 and 1/32) for a = 3 (Fig. 2(b)). 

In the case of q = 1, the recurrence-time^ distribution 
P{T) for smallej: a exhibits a weak peak at T/T --^0.1 and a 
shoulder at T/T ~ 1, with a tail of the distrilDution showing 
an exponential behavior at longer T. In the case of a = 3, 
a smaller-T peak in P{T) becomes more pronounced, while 
a larger-T shoulder gets almost lost. While the computed 
P{T) at larger T does not much depend on the grid spacing 
a, the smaller-T peak continues to shift to smaller T as a 
decreases. We note that, in both cases of a = 1 and 3, the 
basic feature of the distribution is robust against the change 
of the range parameter r. 



3. The non-viscosity model 

In this section, we analyze the continuum limit of the ID 
non-viscosity BK model with rj = 0. We show the results of 
our numerical simulations for various observables, including 
the recurrence-time distribution, the time-correlation func- 
tion of seismic events before and after the mainshock, the 
time development of the spatial correlation function of seis- 
mic events before the mainshock, and the time development 
of the magnitude distribution function before and after the 
mainshock, with varying the block-size parameter a in the 
range 1 > a > 1/32. We set a either a = 1 or a = 3. 
As mentioned, a = 1 corresponds to the "near-critical" case 
while the latter a = 3 corresponds to the "supercritical" 



3.1. The magnitude distribution 

The magnitude distributions Rin) of the ID non- viscosity 
BK model is shown in Figs. 1(a) and (b) for the cases of 
Q = 1 and 3, respectively. The grid spacing a is varied 
in the range 1 > a > 1/32. One can see from the figures 
that the qualitative features of i?(/i), i.e., a "near-critical" 
feature in the case of a = 1 and a "supercritical" feature 
in the case of q — 3, are preserved even in the continuum 
limit. Meanwhile, in both cases of a = 1 and a = 3, the 



3.3. Time correlations of events associated with the 
mainshock 

In Fig. 3, we show the time correlation function between 
large events (mainshock) and events of arbitrary sizes for the 
cases of a = 1 (a) and of a = 3 (b) . In these figures, we plot 
the mean number of events of arbitrary sizes, dominated 
in number by small events, occurring within the distance 
r = 7.5 from the epicenter of the mainshock before (i < 0) 
and after (i > 0) the mainshock, where the occurrence of the 
mainshock is taken to be the origin of the time t — 0. The 
average is taken over all large events of their magnitudes of 
fj. > fic = 2. The number of events are counted here with 
the time bin of AW = 0.02. 

In both cases of a = 1 and 3, qualitative features of the 
time correlation of the original discrete model, including a 
prominent seismic acceleration before the mainshock [Shaw 
et at, 1992; Mori and Kawamura, 2006], persists in the con- 
tinuum limit. Seismic suppression after the mainshock is 
more pronounced in a = 3 than in a = 1. The data of 
Q = 3 do not seem to converge to an asymptotic form even 
at our smallest grid spacing a — 1/32. 

3.4. Spatial correlations of events before the 
mainshock 

In this subsection, we examine the continuum limit of 
the spatial seismic correlations before the mainshock. Our 
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previous studies on the discrete BK model [Mori and Kawa- 
mura, 2005; 2006; 2008a] has revealed the occurrence of the 
doughnut-like quiescence phenomenon prior to the main- 
shock, i.e., though the frequency of small events are gen- 
erally enhanced preceding the mainshock at and around 
the epicenter of the upcoming mainshock, the frequency of 
smaller events is suppressed just before the mainshock in a 
close vicinity of the upcoming mainshock, while it continues 
to be enhanced in the surrounding blocks. This phenomenon 
closely resembles the "Mogi doughnut" [Mogi, 1969; 1979, 
Scholz, 2002]. The spatial range where the quiescence was 
observed was narrow, only of a few blocks. The cause of 
this doughnut-like quiescence was ascribed to the one-block 
events [Mori and Kawamura, 2006; 2008a]. 

Then, a natural question to be addressed is whether the 
doughnut-like quiescence observed in the discrete BK model 
survives the continuum limit, or it is a phenomenon intrin- 
sically originated from the short cut-off length scale of the 
model. 

In order to answer this question, we calculate the spa- 
tial seismic correlation functions before the mainshock, with 
systematically varying the grid spacing a. The seismic cor- 
relation between the mainshock of /i > /i^ = 2 and the 
preceding events of arbitrary size, dominated in number by 
small events, are calculated for several time periods before 
the mainshock. Fig. 4 exhibits the time-dependent spatial 
correlation functions in the case of q = 1, for a larger grid 
spacing a = 1/4 (a) and for a smaller grid spacing a = 1/32 
(b), while Fig. 5 exhibits the same quantities in the case of 
Q = 3. In order to see the dependence on the grid spacing 
a more directly, we show in Fig. 6 the spatial seismic cor- 
relation functions in a shorter time period oi ti^ — ^ 0.01 
for various a, in the cases of a = 1 (a) and of a = 3 (b). 
As is clear from Fig. 6, the spatial range of the quiesence 
gets narrower for smaller grid spacing a. Namely, while the 
doughnut-like quiescence is detectable for smaller grid spac- 
ings of a < 1, as the grid spacing a gets significantly smaller, 
however, the spatial range of the quiescence gets narrower, 
tending to vanish for small enough a. See insets of Fig. 6 
where the peak position of the event frequency representing 
the quiesence scale is plotted as a function of the grid spac- 
ing a. Apparently, the peak position tends to zero in the 
continuum limit a — * 0. This observation strongly suggests 
that the doughnut-like quiescence might vanish altogether 
in the continuum limit a ^ 0. 

3.5. The time-dependent magnitude distribution 
before the mainshock 

In the original discrete BK model, a significant change 
of the magnitude distribution is often observed preceding 
the mainshock [Mori and Kawamura, 2005; 2006; 2008a]. 
In the ID short-range BK model, in particular, a signifi- 
cant suppression of large events was observed just before 
the mainshock, leading to an apparent increase of the effec- 
tive B-value. In this subsection, we examine the continuum 
limit of such temporal evolution of the magnitude distribu- 
tion before the mainshock. 

In Figs. 7 and 8, we show the "time-resolved" local mag- 
nitude distributions for several time periods before the main- 
shock in the cases of a = 1 (Fig. 7) and of q = 3 (Fig. 8). 
Figs, (a) and (b) represent the cases of a larger grid spacing 
a = 1/4 and a smaller grid spacing a = 1/16, respectively. 
Only events with their epicenters lying within the distance 
r = 7.5 from the upcoming mainshock of /i > /ic = 2 are 
counted here. As can be seen from these figures, in both 
cases of a = 1 and 3, the continuum limit seems to afi'ect 
only slightly the temporal evolution of the magnitude dis- 
tribution. Even in the continuum limit, the suppression of 
large events is observed before the mainshock, leading to an 
apparent increase of the effective _B-value. 

3.6. The time-dependent magnitude distribution 
after the mainshock 

Our previous study showed that, in contrast to the behav- 
ior before the mainshock, the magnitude distribution of the 



discrete ID BK model after the mainshock showed very lit- 
tle temporal evolution [Mori and Kawamura, 2006; 2008a]. 
Here, we examine how the continuum limit affects this prop- 
erty after the mainshock. 

In Figs. 9 and 10, we show the time-resolved local mag- 
nitude distributions after the mainshock for the case a — 1 
(Fig. 9) and of a = 3 (Fig. 10). Figs, (a) and (b) repre- 
sent the cases of a larger grid spacing a = 1/4 and a smaller 
grid spacing a — 1/16, respectively. Other conditions are 
the same as those in Figs. 7 and 8. As can be seen from 
these figures, in both cases of a = 1 and 3, the magnitude 
distribution in the continuum limit shows little-to-no tem- 
poral evolution after the mainshock. Aftershock sequence 
obeying the Omori law is never realized in the model even 
in its continuum limit. 

4. The viscosity model 

In this section, we show the results of our numerical sim- 
ulations on the ID BK model with a nonzero viscosity term 
rj — 0.02 for various observables studied in the previous 
section. As in the previous section, the velocity-weakening 
parameter a is taken to be either a = 1 or a = 3, whereas 
the parameter a is fixed to a = 0.01. 

4.1. The magnitude distribution 

The computed magnitude distributions R{fi) are shown in 
Figs. 11(a) and (b) for the cases of a = 1 and 3, respectively. 
As can be seen from the figures, a nonzero viscosity tends 
to weaken the GR character of the magnitude distribution. 
Namely, it weakens the linearity of the R{fi) curve, causing 
an eminent down-bending behavior. The main cause of this 
deviation from the GR law comes from the supression of 
events at smaller magnitudes in the viscousity model rela- 
tive to those in the non- viscousity model. Indeed, the viscos- 
ity tends to make the relative displacement of neighboring 
blocks being smoothe, enhancing the correlated motion of 
neighboring blocks. As a result, the frequency of smaller 
events of one or a few blocks is considerably reduced in the 
viscous model, which causes the observed deviation from the 
GR law at smaller magnitudes. 

Notable difference from the non-viscous case shown in 
Fig. 1 is that R{fi) seems to converge to an asymptotic 
form for smaller a in both cases of q = 1 and 3, except that 
the minimum magnitude continuously gets lower as the grid 
spacing a gets smaller. A similar result was already reported 
by [Shaw, 1994]. 

The small-length cut-ofi' scale e as given by Eq.(9) is es- 
timated here to be e ~ 0.44 and 0.26 for a — 1 and 3, 
respectively. As can be seen from Figs. 11(a) and (b), R{fi) 
converges to an asymptotic form for the a-values smaller 
than a ~ 1/4 and 1/8 for a = 1 and 3, respectively, which 
is consistent with the expected condition of the continuum 
limit a < e. 

In Fig. 12, the viscosity tj dependence of -R(/i) is shown 
with fixing a — 1/4. In both cases of a = 1 and 3, as 
the viscosity gets larger, weights of larger events are signifi- 
cantly suppressed. In the case of a = 1 shown in Fig. 12(a), 
_R(/i) exhibits a down-bending curvature, while in the case 
of Q = 3 shown in Fig. 12(b), the peak of R{i-i) at a larger 
magnitude shifts toward a smaller magnitude. 

4.2. The local recurrence-time distribution 

In Fig. 13, we show on a log-log plot the local recurrence- 
time distribution function P{T) for the cases of a = 1 (a) 
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and a = 3 (b). Here, we consider large events of their mag- 
nitude IJ. > fie ~ 2, and compute the local recurrence time 
T with the range parameter r = 7.5. The recurrence time 
is normalized by its mean T, which is Tv = 2.22, 1.85, 1.88, 
1.82, 1.82, 1.85 and 1.85 (respectively for a = 1, 1/4, 1/8, 
1/12, 1/16, 1/24 and 1/32) for a = 1 (Fig. 13(a)), and 
Tu = 39.8, 1.44, 2.01, 2.03, 2.04, 2.06 and 2.03 (respectively 
for a = 1, 1/4, 1/8, 1/12, 1/16, 1/24 and 1/32) for a = 3 
(Fig. 13(b)). 

The behavior of the recurrence-time distribution P(T) in 
the viscous case is qualitatively similar to that in the non- 
viscous case, although the convergence is much faster here 
in the viscous case. In the case of q — 1, the asymptotic 
P(T) exhibits a weak peak at T/T ~ 0.1 and a shoulder at 
T/T ~ 1, with a tail of the distribution showing an expo- 
nential behavior at longer T. In the case of a = 3, as can be 
seen from Fig. 13(b), the asymptotic P{T) exhibits a peak 
at T ~ 0.2. 

4.3. Time correlations of events associated with the 
mainshock 

In Fig. 14, we show the time correlation function between 
large events (mainshock) and events of arbitrary sizes, dom- 
inated in number by small events, for the cases of q = 1 (a) 
and a = 3 (b). The other conditions are the same as in Fig. 
3 in the non- viscous case. 

Again, with decreasing the grid spacing a, the results 
converge to the continuum limit fairly quickly. The obtained 
behavior of the time correlation looks more or less similar to 
the one in the non- viscous case, although the convergence is 
much faster here in the viscous case even including the case 
of a = 3. 

4.4. Spatial correlations of events before the 
mainshock 

In this subsection, we examine the spatial seismic correla- 
tions before the mainshock. The calculational conditions are 
the same as in Figs. 4-6 in the viscous case. The computed 
spatial seismic correlation function immediately before the 
mainshock in the shorter time period oi tv = ^ 0.01 is 
shown in Fig. 15 for various grid spacings a. The frictional 
parameter is either a = 1 (a) or a = 3 (b). 

The doughnut-like quiescence is appreciable even for 
smaller grid spacings of a < 1 for both cases of q = 1 and 3. 
As in the non-viscous case studied in the previous section, 
as the grid spacing a gets smaller, the spatial range of the 
quiescence gets narrower, tending to vanish for small enough 
a: See the insets of Fig. 15. This observation strongly sug- 
gests again that the doughnut-like quiescence might vanish 
altogether in the continuum limit a — > 0. 

Hence, the doughnut-like quiescence observed in the dis- 
crete BK model is likely to be a phenomenon closely re- 
lated to the short-length cut-off scale of the model. This 
seems fully consistent with our observation in the previous 
papers that the one-block events are responsible for the ob- 
served doughnut-like quiescence [Mori and Kawamura, 2006; 
2008a]. 

4.5. The time-dependent magnitude distribution 
before the mainshock 

In Figs. 16 and 17, we show the "time-resolved" local 
magnitude distributions for several time periods before the 
mainshock. The case of a = 1 is shown in Fig. 16, and the 
case of Q = 3 in Fig. 17. Figs, (a) and (b) represent the 
cases of a larger grid spacing a = 1/4 and a smaller grid 
spacing a — l/l6, respectively. The calculational conditions 
are the same as in Figs. 7 and 8 of the non-viscous case. 

In both cases of a = 1 and 3, as can be seen from Figs. 
16 and 17, the temporal evolution of R{fi) prior to the main- 
shock is kept qualitatively similar even when the grid spac- 
ing a gets smaller. Namely, on approaching the mainshock, 
weights of larger events are significantly suppressed. Similar 



behavior was observed in Figs. 7 and 8 for the non-viscous 
case. 

We have also studied the "time-resolved" local magnitude 
distributions after the mainshock (the results not shown 
here). As in the non- viscous case shown in Figs. 9 and 10, 
the magnitude distribution hardly changes on approaching 
the mainshock even for smaller grid spacing a. 



Overall, the statistical properties of the viscosity BK 
model in the continuum limit turn out to be qualitatively 
similar to those of the non-viscosity BK model in the contin- 
uum limit, although the "critical" aspect, e.g., the GR-like 
behavior, is weakened somewhat. By contrast, the conver- 
gence to the asymptotic continuum limit is much faster in 
the viscous case than in the non-viscous case. The man- 
ner of this convergence is well controlled by the length scale 
given by Eq.(3). 

5. Summary and discussion 

Spatiotcmporal correlations of the one-dimonsional 
spring-block (Burridge-Knopoff) model of earthquakes, ei- 
ther with or without the viscosity term, are studied by 
means of numerical computer simulations. The continuum 
limit of the model is examined by systematically investigat- 
ing the model properties with varying the block-size param- 
eter a toward a ^ 0. 

The two types of models are analyzed, i.e., the ID BK 
model with the Kelvin viscosity and the one without the 
Kelvin viscosity. The naive continuum limit of the BK model 
without the viscosity is known to be problematic in that the 
dynamics becomes pathological. The Kelvin viscosity term 
is introduced so that the model dynamics possesses a sensi- 
ble continuum limit. The added viscosity term introduces a 

small-length cut-off scale e = tt ^ / -2- , which indeed turns out 
to control the continuum limit of the model. It also turns 
out that the viscosity affects the statistical properties of the 
model somewhat, making large earthquakes smaller and en- 
hancing the deviation from the GR law. In other words, the 
viscosity suppresses the critical feature of the model. The 
viscosity tends to make the relative displacement of neigh- 
boring blocks being smoother and enhance the correlated 
motion of blocks. 

Probably reflecting the intrinsic problem associated with 
its pathological dynamics, even the statistical properties of 
the non- viscosity BK model sometimes continues to change 
even at our smallest grid spacing studied. In contrast, those 
of the viscosity BK model are found to converge well to an 
asymptotic continuum limit once the grid spacing is taken 
smaller than the above small-length cut-off scale e. One ob- 
vious consequence of the continuum limit a — > limit is that 
the size of minimum earthquakes gets smaller. It thus ap- 
pears that, in the continuum limit of the model, even an in- 
finitesimal earthquake is possible. By contrast, many of the 
statistical properties of the original discrete BK model are 
kept qualitatively the same even in the continuum limit: For 
example, the model exhibits a seismic acceleration preceding 
the mainshock. The magnitude distribution of the model ex- 
hibits a significant temporal change before the mainshock, 
while it exhibits only a negligible change after the main- 
shock. This observation in turn gives some guarantee that 
the original discrete BK model provides a reasonable de- 
scription of the continuum fault. 

One notable difference between the original discrete 
model and the corresponding model in the continuum limit is 
the existence/non-existence of the doughnut-like quiescence 
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just before the mainshock. Large events of the original dis- 
crete BK model exhibits a seismic acceleration before the 
mainshock, which is accompanied by a doughnut-like qui- 
escence occurring just before the mainshock in close vicin- 
ity of the upcoming mainshock. It is found that in both 
viscous and non-viscous cases, as the grid spacing a gets 
smaller, the spatial range of the doughnut-like quiescence 
becomes narrower, and the doughnut-like quiescence might 
vanish altogether in the continuum limit. The doughnut- 
like quiescence observed in the discrete BK model is then a 
phenomenon closely related to the short cut-off length scale. 

This observation might have some implications to real 
seismicity. While the real crust is obviously a continuum, 
it is not so obvious whether there exists an intrinsic short- 
length cut-off scale in real seismicity or not. In any case, in 
real earthquakes, the "Mogi-doughnut" is occasionally re- 
ported to occur [Mogi, 1969; 1979, Scholz, 2002], ahhough 
to elucidate its statistical significance is often not easy. Our 
present result might suggest that, if the real crust possesses a 
cut-off length scale, the "Mogi- doughnut" quiescence might 
occur at such a length scale. In other words, spatial inho- 
mogeneity might be an essential ingredient for the "Mogi- 
doughnut" to occur. 

For the future, it is desirable to extend the present analy- 
sis of the continuum limit of the ID BK model to the corre- 
sponding 2D model. Since taking the a — > limit necessarily 
means taking the large-size limit, performing such analysis 
in 2D is computationally demanding. It is also desirable 
to perform a similar analysis for the BK model based on 
the rate- and state-dependent friction law [Morimoto and 
Kawamura, 2007], which is the friction law employed in the 
study of Rice [Rice, 1993]. This friction law naturally pro- 
vides a length scale into the problem, even without invok- 
ing the Kelvin viscosity. Thus, in order to fully resolve the 
question raised by Rice, it would be necessary to perform a 
similar analysis for the BK model with the rate- and state- 
dependent friction law. We leave such interesting extension 
of the present analysis to a future task. 
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Figure 1. The magnitude distribution -R(/i) of earth- 
quake events of the ID non- viscosity BK model (77 = 0). 
The dimensionless grid spacing a is varied in the range 
1 > a > 1/32. Figs, (a) and (b) represent the cases of 
a — 1 and 3, respectively. The parameter a is fixed to 
a — 0.01. The system size is L = aN = 200, except 
the case of a = 1/32 of a = 1 where L is taken to be 
L — aN = 400 to circumvent the more severe finite-size 
effect. 
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Figure 2. The local recurrence-time distribution func- 
tion P{T) of the ID non- viscosity BK model (t; = 0), with 
varying the dimensionless grid spacing a. Figs, (a) and 
(b) represent the cases of a = 1 and 3, respectively. The 
parameter a is fixed to ct = 0.01. The main panels repre- 
sent the log- log plots of P{T), while the insets represent 
the semi-logarithmic plots including the tail part of the 
distribution. The mean recurrence time T is Tv = 2.22, 
0.56, 2.30, 2.67, 2.59 and 3.20 (respectively for a = 1, 1/4, 
1/8, 1/12, 1/16 and 1/24) for a = 1, and Tv = 41.4, 1.72, 
2.76, 3.90, 4.81, 6.09 and 6.72 (respectively for a = 1, 1/4, 
1/8, 1/12, 1/16, 1/24 and 1/32) for a = 3. 
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Figure 3. The time correlation function of the ID 
non-viscosity BK model {rj = 0) between large events of 
^c ~ 2 (mainshock) occurring at time t = and events 
of arbitrary sizes (dominated in number by small events) 
occurring at time t. The dimensionless grid spacing a is 
varied in the range 1/4 > a > 1/32. Fig. (a) represents 
the case of a = 1, while Fig. (b) represents the case of 
a — 3. The parameter a is fixed to a — 0.01. Events 
of arbitrary sizes occurring within the distance r — 7.5 
from the epicenter of the mainshock are counted. The 
negative time t < represents the time before the main- 
shock, while the positive time f > represents the time 
after the mainshock. The average is taken over all large 
events of its magnitude /i > /Xc = 2. The system size is 
L = aN = 200. 
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Figure 4. Event frequency for several time periods be- 
fore the mainshock oi fi > fj,c = 2 oi the ID non-viscosity 
BK model (77 — 0) plotted versus r, the distance from the 
epicenter of the upcoming mainshock. The frictional pa- 
rameters are a = 1 and a = 0.01. The dimensionless grid 
spacing a is a = 1/4 (a), and a — 1/32 (b). The system 
size is L = aN — 200. The insets represent similar plots 
at longer times. 
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Figure 5. Event frequency for several time periods be- 
fore the mainshock of /i > /ic = 2 of the ID non-viscosity 
BK model (r; = 0) plotted versus r, the distance from the 
epicenter of the upcoming mainshock. The frictional pa- 
rameters are a = 3 and a = 0.01. The dimensionless grid 
spacing o is o = 1/4 (a), and a — 1/32 (b). The system 
size is L = aN = 200. The insets represent similar plots 
at longer times. 
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Figure 6. Event frequency in the time period tv = Q ^ 
0.01 immediately before the mainshock of ^ > /Xc = 2 
of the ID non-viscosity BK model (r\ = 0) plotted ver- 
sus r, the distance from the epicenter of the upcoming 
mainshock. Fig. (a) represents the case of a = 1, while 
Fig. (b) represents the case of a = 3. The dimensionless 
grid spacing a is varied in the range 1/32 < a < 1/4. 
The parameter a is fixed to a = 0.01 The system size 
is 1/ = aN = 200. The insets represent the peak posi- 
tion of the event frequency, corresponding to the range 
of the doughnut-Uke quiescence, as a function of the di- 
mensionless grid spacing a. In both cases of a = 1 and 
3, the doughnut-hke quiescence appears to vanish in the 
continuum hmit a ^ 0. 
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Figure 7. The local magnitude distribution for several 
time periods before the mainshock of /x > /Xc = 2 of the 
ID non- viscosity BK model (»7 = 0). The frictional pa- 
rameters are a = 1 and a = 0.01. Figs, (a) and (b) 
represent the cases of a = 1/4 (a) and a — 1/16 (b), 
respectively. Events whose epicenter lies within the dis- 
tance r = 7.5 from the epicenter of the upcoming main- 
shock are counted. The system size is L = aN = 200. 
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Figure 8. The local magnitude distribution for several 
time periods before the mainshock of /x > /Xc = 2 of the 
ID non- viscosity BK model (77 = 0). The frictional pa- 
rameters are a = 3 and a = 0.01. Figs, (a) and (b) 
represent the cases of a = 1/4 (a) and a = 1/16 (b), re- 
spectively. Events whose epicenter lie within the distance 
r = 7.5 from the epicenter of the upcoming mainshock 
are counted. The system size is L = aN — 200. 
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Figure 9. The local magnitude distribution for several 
time periods after the mainshock of /i > /Xc = 2 of the 
ID non- viscosity BK model (r; = 0). The frictional pa- 
rameters are a = 1 and a — O.OL Figs, (a) and (b) 
represent the cases of a = 1/4 (a) and a = 1/16 (b), re- 
spectively. Events whose epicenter lie within the distance 
r — 7.5 from the epicenter of the preceding mainshock are 
counted. The system size is L = aN = 200. 
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Figure 10. The local magnitude distribution for sev- 
eral time periods after the mainshock oi fj. > fic ~ 2 oi 
the ID non- viscosity BK model {rj = 0). The frictional 
parameters are a = 3 and a — 0.01. Figs, (a) and (b) 
represent the cases of a = 1/4 (a) and a — 1/16 (b), re- 
spectively. Events whose epicenter lie within the distance 
r — 7.5 from the epicenter of the preceding mainshock are 
counted. The system size is L = aN = 200. 
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Figure 11. The magnitude distribution R{n) of earth- 
quake events of the ID viscosity BK model {rj — 0.02) 
with a — 0.01. The dimensionless grid spacing a is var- 
ied in the range 1 > a > 1/32. Figs, (a) and (b) represent 
the cases of a = 1 and 3, respectively. The system size is 
L = aN ^ 200. 
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Figure 12. The magnitude distribution R{n) of earth- 
quake events of the ID viscosity BK model for various 
values of the viscosity parameter rj. Figs, (a) and (b) 
represent the cases of a = 1 and 3, respectively. The 
other parameters are fixed to a = 1/4 and a = 0.01. The 
system size is L = aN = 200. 
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Figure 13. The local recurrence-time distribution func- 
tion P{T) of the ID viscosity BK model (77 = 0.02), with 
varying the dimensionless grid spacing a. Figs, (a) and 
(b) represent the case of a = 1 and 3, respectively. The 
parameter a is fixed to a — 0.01. The main panels repre- 
sent the log- log plots of P{T), while the insets represent 
the semi-logarithmic plots including the tail part of the 
distribution. The mean recurrence time T is Tv = 2.22, 
1.85, 1.88, 1.82, 1.82, 1.85 and 1.85 (respectively for 
a = ]_, 1/4, 1/8, 1/12, 1/16, 1/24 and 1/32) for a = 1, 
and Tv = 39.8, 1.44, 2.01, 2.03, 2.04, 2.06 and 2.03 (re- 
spectively for a = 1, 1/4, 1/8, 1/12, 1/16, 1/24 and 1/32) 
for a = 3. 
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Figure 14. The time correlation function of the ID 
viscosity BK model {i] — 0.02) between large events of 
fic ~ 2 (mainshock) occurring at time f = and events 
of arbitrary sizes (dominated in number by small events) 
occurring at time t. The dimensionless grid spacing a is 
varied in the range 1/4 > a > 1/32. Fig. (a) represents 
the case of a = 1, while Fig. (b) represents the case of 
a — 3. The parameter a is fixed to a = 0.01. Events 
of arbitrary sizes occurring within the distance r = 7.5 
from the epicenter of the mainshock are counted. The 
negative time t < represents the time before the main- 
shock, while the positive time i > represents the time 
after the mainshock. The average is taken over all large 
events of its magnitude /i > /Xc = 2. The system size is 
L = aN = 200. 
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Figure 15. Event frequency in tlie time period ti/ = ^ 
0.01 immediately before the mainshock of /.i > /Xc = 2 of 
the ID viscosity BK model (77 = 0.02) plotted versus r, 
the distance from the epicenter of the upcoming main- 
shock. Fig. (a) represents the case of a = 1, while Fig. 
(b) represents the case of q = 3. The dimensionless grid 
spacing a is varied in the range 1/4 > a > 1/32. The 
parameter a is fixed to cr = 0.01. The system size is 
L = aN = 200. The insets represent the peak position 
of the event frequency, corresponding to the range of the 
doughnut-like quiescence, as a function of the dimension- 
less grid spacing a. In both cases of a = 1 and 3, the 
doughnut-like quiescence appears to vanish in the con- 
tinuum limit a — > 0. 
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Figure 16. The local magnitude distribution for sev- 
eral time periods before the mainshock oi jj, > fj,c = 2 
of the ID viscosity BK model {-q = 0.02). The frictional 
parameters are a = 1 and a = 0.01. Figs, (a) and (b) 
represent the cases of a = 1/4 (a) and a — 1/16 (b), 
respectively. Events whose epicenter lies within the dis- 
tance r = 7.5 from the epicenter of the upcoming main- 
shock are counted. The system size is L = aN = 200. 
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Figure 17. The local magnitude distribution for sev- 
eral time periods before the mainshock of /x > pc = 2 
of the ID viscosity BK model {r/ = 0.02). The frictional 
parameters are a = 3 and a — 0.01. Figs, (a) and (b) 
represent the cases of a = 1/4 (a) and a — 1/16 (b), 
respectively. Events whose epicenter lies within the dis- 
tance r = 7.5 from the epicenter of the upcoming main- 
shock are counted. The system size is L = aN = 200. 



